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Abstract 

Using the Wang-Landau Monte Carlo method, we study the antiferromagnetic (AF) three-state 
Potts model with a staggered polarization field on the square lattice. We obtain two phase tran- 
sitions; one belongs to the ferromagnetic three-state Potts universality class, and the other to the 
Ising universality class. The phase diagram obtained is quantitatively consistent with the transfer 
matrix calculation. The Ising transition in the large nearest-neighbor interaction limit has been 
made clear by the detailed analysis of the energy density of states. 
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I. INTRODUCTION 



The g-state Potts model is one of the basic models for studying phase transitions [1]. 
The properties of the antiferromagnetic (AF) Potts models are more complex than those of 
the ferromagnetic (F) ones. The phase transitions of the AF Potts models depend heavily 
on the number of states q, the details of lattice structure, etc. The AF three-state (q = 3) 
Potts model on a square lattice exhibits a second-order transition at T=0 with the Gaussian 
criticality [2, 3]. While this model with only the nearest-neighbor (NN) interactions has 
been studied in detail [2-9] , the effect of the next-nearest-neighbor (NNN) interactions has 
many interesting problems to explore [4, 10]. Quite recently, the present authors have 
studied the square-lattice AF three-state Potts model with a staggered polarization field 
[10]. By the use of the exact diagonalization calculation of the transfer matrix and the 
phenomenological renormalizat ion-group analysis, two types of phase transitions have been 
discussed in connection with a field theoretical argument. The crossover behavior from the 
AF to the F three-state Potts criticality, which was proposed by Delfino [11], has been 
confirmed. 

The exact diagonalization calculation and the Monte Carlo simulation play a complemen- 
tary role in the numerical study. The energy levels obtained by the exact diagonalization are 
highly accurate, but the tractable size is limited. On the other hand, although the statisti- 
cal errors are unavoidable because of the sampling process, the Monte Carlo simulation can 
deal with larger systems. Moreover, the latter can easily study the spin configuration, the 
probability distribution of the order parameter, etc. Recently, several attempts have been 
proposed for the Monte Carlo algorithms to directly calculate the energy density of states 
(DOS), such as the multicanonical method [12, 13], and the Wang-Landau method [14]. 

In this paper, we study the square-lattice AF three-state Potts model with a staggered 
polarization field by using the Wang-Landau Monte Carlo method. We calculate the energy 
DOS precisely, and study the phase transitions of the model. Using the finite-size scaling 
(FSS), we investigate the critical properties. In Sec. II, the model and the calculation method 
are described. The results for phase transitions are given in Sec. III. The final section is 
devoted to a summary and discussions. 
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II. MODEL AND CALCULATION METHOD 



We treat the AF three-state Potts model with the staggered polarization field on the 
square lattice A, whose Hamiltonian is given as 

h = *E^-^£(-i)'\W (i) 

(j,k) [j,k] 

Here, <jj = 0, 1, 2, and J±, J<i > 0. The first sum is performed over the whole NN pairs (j, k), 
and the second sum over the whole NNN pairs [j,k]. Here, (— 1) J = ±1 for j in the even 
(odd) sublattice A±. The second term is the staggered polarization field term, which breaks 
the sublattice symmetry. The NN and NNN interactions are schematically shown in Fig. 1 
for convenience. We also give an example of the ground-state configuration there. 

In order to obtain precise numerical information we use the Wang-Landau Monte Carlo 
algorithm [14] because of the following reasons. This type of extended ensemble methods 
do not suffer from the problem of the critical slowing down near the second-order phase 
transitions. The cluster algorithms are not so efficient for the systems with NNN couplings. 
We use the information of the energy DOS for the detailed study of the phase transitions. In 
the Wang-Landau method, a random walk in energy space is performed with a probability 
proportional to the reciprocal of the DOS, l/g(E), which results in a flat histogram of energy 
distribution. Since the DOS is not known a priori, it is iteratively updated as 

In <7(£)- In </(£?)+ In/, (2) 

every time a random walker visits a state with energy E. A large modification factor / is 
introduced to accelerate the diffusion of the random walk in the early stage of the simulation, 
and it is gradually reduced to unity by checking the 'flatness' of the energy histogram. We 
set the final value of (In/) as 1CT 8 following the original paper by Wang and Landau [14], 
and then measure the energy dependence of a quantity Q, that is, Q(E). In measuring 
Q{E), the DOS g(E) is fixed as the final one. Finally we calculate the thermal average (Q)t 
at the temperature T using the energy DOS and the Boltzmann weight as 

fns E E Q(E)g(E)e~ E / T 

{Q)T - Y. E 9{E)e-W ' (3) 
where the Boltzmann constant has been included in the definition of T. 

We make simulations for several sets of J2/J1, we treat the system sizes N = L x L up 
to L = 64. For the energy range of the random walk we do not cover the whole possible 
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energy space to save computation time. The lowest energy is taken as the ground-state 
energy, E = —J2N, but we set the highest energy as that takes the maximum DOS, in 
other words, the energy at the infinite temperature. This energy is (2/3)J 1 N. We make 
the measurement for 4 x 10 6 Monte Carlo steps per spin after the final energy DOS is 
obtained using the Wang-Landau process. We perform 64 independent runs for each system 
size in order to get better statistics and to estimate statistical errors. We estimate the 
statistical errors from 64 independent calculations. They could be underestimated if there 
are systematic errors due to the Wang-Landau method. 

III. RESULTS 

A. Specific heat 

First, we present the data for the specific heat of the AF three-state Potts model with 
the staggered polarization field on the square lattice. We plot the temperature dependence 
of the specific heat for J 2 /Ji = 1/2, 1, and 2 in Fig. 2. The temperature is denoted in units 
of T/Ji from now on unless specified else. The system sizes are L=16, 24, 32, 48, and 64. 
The statistical errors are within the width of lines. 

We see two peaks in the specific heat for each J 2 /Ji, which indicates the existence of two 
phase transitions clearly. The peak value for the high-temperature phase transition increases 
rapidly with the system size, which suggests a positive specific-heat exponent a. The peak 
value for the low-temperature transition, on the contrary, increases gradually with the size. 
It is difficult, however, to determine whether this increase is a logarithmic divergence or a 
power- law one with very small ajv for these system sizes. 

B. Order parameter 

To investigate the behavior of the phase transitions in more detail, we consider the order 
parameters. For high-temperature phase transition, we look at the staggered magnetization 
for the AF Potts model, which is given by 

M 2 = S 2 + S, 2 + S 2 2 , (4) 
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where 

^=^£(-1)^, € = 0,1,2. (5) 

JV jeA 

The order parameter M s takes non-zero values if the S3 symmetry associated with the global 
permutations of three Potts states is broken. At the ground state of the AF three-state Potts 
model with the staggered polarization field, which is shown in Fig. 1, this order parameter 
takes the value of M s 2 = 3/8; we should note that the ground state is sixfold degenerate. 
The temperature dependence of (M s 2 ) is plotted in Fig. 3(a). The value of J2/J1 is 1/2, 1, 
and 2; for each J2/J1, the data for L=16, 24, 32, 48, and 64 are shown. We see that M s 
grows below the temperature that gives a peak of the specific heat, which indicates that M s 
is an appropriate order parameter for the high-temperature phase transition. The F order 
develops on the sublattice A + . 

For the quantitative analysis of the phase transition, we use the moment ratio, 
(M S 4 )/(M S 2 ) 2 , which is essentially the same as the Binder ratio [15]. We plot the mo- 
ment ratio for the high-temperature phase transition in Fig. 3(b). Curves with different 
sizes cross at a single point if the corrections to FSS are negligible. The crossings of our 
data are very good. From the crossing point, we estimate the critical temperature. We can 
also estimate the critical exponent v by the FSS analysis, 

(M S 4 )/(M S 2 ) 2 = fitL^), (6) 

where t = (T — T c )/T c . We show the FSS plot of the moment ratio for J 2 /Ji=l, for example, 
in Fig. 4(a). We see a very good FSS. We can also estimate the critical exponent f3/v using 
the FSS relation, 

(M s 2 ) T=Tc ~ L-W\ (7) 

We tabulate the estimates of T c , v, and f3/v for the high-temperature phase transition 
in Table I. The results for J 2 /Ji = 1/4 and 4 are also given there. The numbers in the 
parentheses denote the uncertainty in the last digits. We use the least-square fitting for 
the FSS estimates without considering the corrections to FSS. The statistical errors are 
estimated from 64 independent runs. We see from Table I that v is consistent with the F 
three-state Potts value of v — 5/6 = 0.833. We also find that f3/v is compatible with the 
F three-state Potts value of (3/v = 2/15 = 0.133. The estimates of the moment ratio at 
T = T c are also given in Table I. This value is consistent with that for the F three-state 
Potts model, 1.16 ±0.01 [16]. 
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Here, the corrections to FSS have not been considered explicitly for the estimate of the 
critical temperature and critical exponents. The J-ij J\ dependences of the estimated critical 
exponents and the moment ratio at T = T c are small, but the deviations from the values for 
the F three-state Potts model become larger for small J-ij J\. This behavior is prominent 
for the moment ratio at T = T c , which is more accurate than the critical exponents. This is 
consistent with the fact that the system on the sublattice A + is nothing but the F three-state 
Potts model for large J2/J1 limit, which determines the renormalizat ion-group flow [10]. 

Next consider the low-temperature phase transition. The S 3 symmetry is already broken 
in the sublattice A + ; we should consider the symmetry breaking within the sublattice A_. 
Then, we may consider the following quantity as the order parameter for the low-temperature 
phase transition: 

m s 2 = s 2 + Sl 2 + s 2 2 , (8) 

where 

st = ^ EH)*^- € = 0,1,2. (9) 

We only look at the spins on the sublattice A_. The sublattice A_ forms a \pl x \[2 square 
lattice, and we divide the sublattice A_ into two further sublattices. For these further 
sublattices, we take (— l) k = ±1 depending on the values of k, even or odd. This order 
parameter essentially represents the AF Ising order in the sublattice A_. At the ground 
state the low-temperature order parameter, Eq. (8), takes the value of m s 2 = 1/8. 

In Fig. 5(a) we plot the order parameter of the low-temperature phase transition, Eq. (8). 
This time, m s grows below the temperature that gives a lower peak of the specific heat; we 
see that m s is an appropriate order parameter for the low-temperature phase transition. We 
can consider the moment ratio associated with the low-temperature order parameter; M s 
will be replaced by m s in Eq. (6). This moment ratio is plotted in Fig. 5(b). Again, we see 
the crossing of the curves with different sizes. Using the FSS analysis, we estimate T c , v, 
and (3/ v, and they are tabulated in Table I. The FSS plot of the moment ratio for J 2 / J±=l 
is given in Fig. 4(b), for example. The estimated exponents shown in Table I suggest that 
the exponents v and /3/v are consistent with the two-dimensional (2D) Ising values, 1 and 
1/8=0.125, respectively. The moment ratio at T = T c is also given in Table I. The Ising 
value for the finite system with aspect ratio 1:1 is rigorously calculated as 1.1679229(47) 
[17]. The coincidence with this value is very good. 
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The J2/J1 dependences of the estimated critical exponents and the moment ratio at 
T = T c are very small, but the deviations from the Ising values become slightly larger for 
large J 2 / Ji, which is the opposite direction from the case of the high-temperature transition. 

C. Phase diagram 

From the list of the estimates of the two critical temperatures given in Table I for J 2 j J\ = 
1/4, 1/2, 1, 2, and 4, we discuss the phase diagram. In Fig. 6 we plot the phase diagram in the 
parameter space of (u,v) = (e~ Jl ^ T , 1 — e' j2 ^ T ), which is the same as the previous transfer- 
matrix study [10]. The trajectories of J 2 /Ji = l/4, 1/2, 1, 2, and 4 are shown by dotted 
curves starting from T = [(u,v) = (0, 1)] to T = 00 [(u,v) = (1,0)]. The filled circle (red) 
and filled square (blue) represent the estimate of high-temperature and low-temperature T c 's 
in the present study, respectively, and they are compared with the previous estimates by 
the transfer-matrix calculation [10], which are given by open marks. These two results are 
consistent with each other, which shows the reliability of both calculations. In the transfer- 
matrix calculation [10], the characterization of the excitation levels was important. In the 
present MC simulation, we have employed the method to calculate the DOS accurately, and 
have made the appropriate choice of the order parameters. With these careful treatments, 
we have obtained the accurate enough phase diagram. 

From the behavior of the corrections, we can see the renormalization-group flow of the 
phase boundary, which is consistent with Zamolodchikov's c-theorem [18]. For the high- 
temperature phase transition, the flow starts from the Gaussian point (u, v) = (0, 0) to the 
F three-state Potts point (u,v) = (1, (3 — v / 3)/2 = 0.6340), which is shown by the double 
circle. On the other hand, for the low-temperature phase transition, the flow starts another 
Gaussian point (u,v) = (1,1) to a point with 2D-Ising criticality on the t>-axis, which is 
shown by the arrow in Fig. 6. 

This Ising transition in the large J\ limit (on the v -axis) is a subtle problem. In Ref . [10] , 
the inverse critical temperature was estimated as Ji/T c = 0.8820, which is slightly larger 
than the Ising value of ln(>/2 + 1)=0.8814 [v = 1 - e' J ^ T =2 - ^=0.5858]. 
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D. Large NN interaction limit 

Here we examine the Ising transition in the large J\ limit. If we consider only the NN 
interaction term in Eq. (1), the energies of the ground states and the first excited states 
are and Ji, respectively. The NNN interaction term takes the energy between —J 2 N 
to J 2 N . Thus, for the case of J 2 j J\ < 1/(2 JV), we have to consider only the ground- 
state configuration for NN interactions. We make a Wang-Landau type MC simulation for 
our model with this restriction. That is, only the ground-state configurations for the NN 
AF three-state Potts model are allowed, and we calculate the energy DOS for the NNN 
interactions. 

In Fig. 7(a), we plot the DOS for this restricted model in the subspace of the ground states 
of the NN AF three-state Potts model, g a {E). If all the spins on the sublattice A + take one 
of the three states (the complete F three-state Potts order) and the spins on the sublattice 
A_ take either one of the other two states, these spin configurations are parts of the ground 
states of the NN AF three-state Potts model, and are called the broken-sublattice-symmetry 
states in the study of the AF Potts models [19]. Then, the spins on the \[2 x y/2 sublattice 
A_ are regarded as the AF Ising model with the NN interactions. The energy DOS of 
this pure Ising model, gb(E), is also shown in Fig. 7(a). The extra three-times degeneracy 
due to three complete F states in the sublattice A + was taken into account. The phase 
transition is determined by the structure of the energy DOS near the Ising critical energy, 
(E/J 2 )/N = -(2 + v / 2)/4 = -0.8536, which is shown by an arrow in Fig. 7. In this region 
two DOS's look close when the logarithmic scale is used. However, we understand the very 
small difference in the critical temperature from the structure of DOS. We plot the ratio of 
g b (E)/g a (E) in Fig. 7(b). Then, we find that this ratio becomes smaller for larger system 
size L. This means that we cannot ignore contributions from the spin configurations which 
are not the pure Ising one. Therefore, we may conclude that although this transition belongs 
to the Ising universality class, the critical temperature, which is not a universal quantity, is 
slightly modified from that of the pure Ising model. Since the number of configurations are 
larger than that for the pure Ising model, the critical temperature T c / J 2 becomes slightly 
lower than that for the pure Ising model. In Fig. 8, we compare the specific heat calculated 
from two DOS's, which give almost the same behavior but there is still a small difference. 

In the MC method to calculate DOS we obtain the relative ratio of DOS's at different 
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energies E\ and E 2 , g(E 1 )/g(E 2 ). The absolute value of g(E) can be obtained with other 
conditions. For the Ising model, as an example, the equation to give the total number of 
states, Yje9{E) = 2 N , is such a condition. There is also a boundary condition if the ground- 
state degeneracy is known. The ground state of the present model is sixfold degenerate; thus 
g(E = —J 2 N)=6. Then, with this condition we can calculate the total number of states, 
J2e9(E). This value is nothing but the ground-state degeneracy for the NN AF three- 
state Potts model, because we restrict ourselves to the subspace of the ground states of the 
AF three-state Potts model. To confirm the reliability of our calculation, we calculate the 
normalized ground-state entropy of the AF three-state Potts model, S/N = [\to.J2e 9{E)]/N, 
as a function of L with this procedure, S(l)/N, and they are tabulated in Table II. More 
direct way of obtaining the ground-state entropy is to calculate the ground-state DOS of the 
three-state Potts model with the condition of J2e 9{E) — 3 N ■ This procedure was employed 
in the study of the three-dimensional AF Potts models [20]. The result of this direct way, 
S(2)/N, is also tabulated in Table II. These two results coincide with each other completely. 
This indicates the effectiveness of the treatment of this subsection, that is, we have applied 
the Wang-Landau method to the present model with the restricted configurations. The 
normalized ground-state entropy for the AF three-state Potts model on the square lattice is 
exactly known as ln(4/3) 3 /2 = 0.431523 for the infinite size limit [21]. The extrapolation of 
the data of S/N shown in Table II as a polynomial of 1/N yields 0.43154 + 1.10 x (1/N), 
which is consistent with the exact value within the statistical errors. We should mention 
that the ground-state entropy of the AF Potts model was extensively studied numerically 
by Shrock and Tsai [22]. 

IV. SUMMARY AND DISCUSSIONS 

We have studied the square-lattice AF three-state Potts model with a staggered polar- 
ization field using the Wang-Landau MC method. We have obtained two phase transitions, 
which belong to the ferromagnetic three-state Potts and Ising universality classes. We 
have confirmed the quantitative consistency of the phase diagram with the transfer-matrix 
study [10]. This consistency shows the reliability of both calculations, the previous transfer- 
matrix calculation and the present MC study. A special attention has been paid to the Ising 
transition in the large J\ limit. The origin of the slight difference of the transition point 
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from the value of the pure Ising model has been made clear by the detailed analysis of the 
energy DOS. As a check of the method, we have calculated the ground-state residual entropy 
for the AF three-state Potts model by two ways. 

The specific heat has singularities at the critical points as shown in Fig. 2. The specific- 
heat amplitudes have information on the crossover behavior [23] . We can study the Gaussian 
to F three-state Potts and the Gaussian to Ising crossovers by the J-ij J\ dependence of the 
specific-heat amplitudes. The detailed study on this crossover will be left to a separate 
study. 

We here make comments on the Wang-Landau method. In calculating the thermal average 
of M s and m s , we have used the microcanonical average M S (E) and m s (E) as in Eq. (3). If the 
energy DOS g(E) is not converged enough, it may cause a systematic error in calculating 
the thermal average. The use of the joint DOS g(E,Q) may help the convergence. We 
have checked for smaller system sizes that the present calculations with the microcanonical 
average and those using the joint DOS give the same results within the statistical errors. 
The second comment is as follows. We may employ a random walk in the space of two 
parameters, J2ij,k) ^aj,a k and J2\j,k](~ l)" ? ^-,o- fc , instead of the single parameter of the total 
energy E. Then, we can get all the information for different J2/J1 from the result of a single 
simulation. Actually, for smaller system sizes (L < 16), the two-parameter random walk 
method works well. However, simulations with fixed J\ and J 2 are more effective for larger 
sizes. 

In this paper, we have considered the staggered polarization field for the NNN interaction. 
The effect of the F NNN interaction is also interesting. In this case, the system has the 
same universality class as the six-state clock model, which yields two Berezinskii-Kosterlitz- 
Thouless transitions [24, 25] . The precise calculation of this problem has been quite recently 
performed by using the level-spectroscopy method [26] based on the exact diagonalization 
of the transfer matrix [27] . A complementary study of the AF Potts model with the F NNN 
interaction using the Monte Carlo method is highly needed, and this research is now in 
progress. 
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FIG. 1: Schematic illustration of the AF three-state Potts model with the staggered polarized 
field. The sites belonging to the sublattice A + (A_) are denoted by red thick (black thin) circles. 
The NNN interactions on the sublattice A + (A_) are F (AF). An example of the ground-state 
configuration is also given. 
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FIG. 2: The specific heat of the AF three-state Potts model with the staggered polarization field 
for J%l J\ = 1/2, 1, and 2. The temperature is plotted in units of T j J\. The system sizes are L= 
16, 24, 32, 48, and 64. The statistical errors are within the width of lines. 
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FIG. 3: The high-temperature order parameter (a) and the ratio of the its moments (b) of the 
AF three-state Potts model with the staggered polarization field for J2/J1 = 1/2, 1, and 2. The 
temperature is plotted in units of TjJ\. The system sizes are L= 16, 24, 32, 48, and 64. The 
statistical errors are within the width of lines. 
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FIG. 4: Scaling plot of the moment ratios of the AF three-state Potts model with the staggered 
polarization field for J2/J1 = 1- We plot both the high-temperature order parameter (a) and the 
low-temperature one (b). The temperature is plotted in units of T/J\. The system sizes are L= 
16, 24, 32, 48, and 64. 
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FIG. 5: The low-temperature order parameter (a) and the ratio of the its moments (b) of the 
AF three-state Potts model with the staggered polarization field for J2/J1 = 1/2, 1, and 2. The 
temperature is plotted in units of TjJ\. The system sizes are L= 16, 24, 32, 48, and 64. The 
statistical errors are within the width of lines. 
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FIG. 6: Phase diagram of the AF three-state Potts model with the staggered polarization field. 
The filled circle (red) and filled square (blue) represent the estimate of high-temperature and low- 
temperature T c 's in the present Monte Carlo study, whereas the open marks represent those by the 
transfer-matrix calculation [10]. The F three-state Potts point (u,v) = (1,(3 — v3)/2 = 0.6340) 
is shown by the double circle, and a point with 2D-Ising criticality on the v-axis is shown by the 
arrow. 
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(E/J 2 )/N 

FIG. 7: (a) The energy DOS of the AF three-state Potts model with the staggered polarization 
field in the subspace of the ground states of the NN AF three-state Potts model, g a (E), and that 
of the pure Ising model, gb(E). The critical energy of the Ising transition point, (E/J2)/N = 
-(2 + \/2)/4 = -0.8536, is shown by an arrow, (b) The ratio of two DOS's, gi ) (E)/g a (E). 
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1 2 

T/J 2 

FIG. 8: Comparison of the specific heat between that for the present model in the subspace of the 
ground states of the NN AF three-state Potts model (solid line) and that for the pure Ising model 
(dotted line). 
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TABLE I: The list of T c , u, (3/v, and moment ratio at T = T c for two phase transitions. The 
temperature is represented in units of T j J\. 



>hhh 


T c 


V 


P/v 


ratio 


high temp. 










1/4 


0.6720(2) 


0.83(1) 


0.105(4) 


1.120(6) 


1/2 


0.9000(3) 


0.83(1) 


0.114(4) 


1.138(6) 


1 


1.3150(5) 


0.83(1) 


0.122(4) 


1.148(6) 


2 


2.1911(5) 


0.83(1) 


0.126(4) 


1.154(6) 


4 


4.0840(5) 


0.82(1) 


0.127(4) 


1.154(6) 


iow temp. 










1/4 


0.2830(2) 


1.00(1) 


0.120(4) 


1.162(4) 


1/2 


0.5642(3) 


1.01(1) 


0.123(4) 


1.165(4) 


1 


1.0657(5) 


1.01(1) 


0.123(4) 


1.164(4) 


2 


1.8228(5) 


0.99(1) 


0.125(4) 


1.164(4) 


4 


3.0240(5) 


0.98(1) 


0.123(4) 


1.155(4) 


TABLE II: The estimate of the residual entropy per spin of the AF three-state Potts model on the 


square lattice. Two procedures are employed to calculate S/N. 


(1) The total number of states for 


the subspace of the £ 


ground states of the AF three-state Potts model is calculated. (2) The direct 


way to calculate the 


ground-state DOS of the three-state Potts 


model is used with the condition 


of Ee9(E)=3 n . 










L 




S(l)/N 




S(2)/N 


16 




0.43569(8) 




0.43575(11) 


24 




0.43345(6) 




0.43353(7) 


32 




0.43259(4) 




0.43262(9) 


48 




0.43198(3) 




0.43202(15) 


64 




0.43176(5) 




0.43185(10) 
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